load packages and functions

source('mytheme.R')
# model version
modelversion = 'gap_log_rstan'
rstanmodelPath = 'modelrlt'
modelPath = paste0(rstanmodelPath, '/models/', modelversion)

Merge the model Result data

1 load data

AllExpData = read.csv(paste0("../data/AllValidData.csv"))
dur <- sort(unique(AllExpData$curDur))

AllExpData$WMSize <- factor(AllExpData$WMSize, labels = c("low", "medium",  "high")) 
# 1: 500ms,  2: 2500, 3 : 2000ms the mean reaction time of WM test 
AllExpData$gap <- factor(AllExpData$gap,  labels = c('short','long', 'not sure'))
AllExpData[which(AllExpData$Exp == 'Exp4'),"Exp"] = "Exp4a"
AllExpData[which(AllExpData$Exp == 'Exp5'),"Exp"] = "Exp4b"

2 Corrct rate

WMCrr2 <- ggplot(meanForPlot, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,
                                  group =Exp, color = Exp, fill = Exp))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  coord_cartesian(ylim = c(0.5, 1)) +
  colorSet5+
  labs(x = "Memory load", y = "Mean accuracy in WM task") +
  theme_new 
WMCrr2

ggsave(paste0(getwd(), "/figures/WMCrr2.png"), WMCrr2, width = 4, height = 4)
### generate WM correct rates
AllExpData$WMCrr <- AllExpData$TPresent == AllExpData$WMRP
m_wmp<- dplyr::group_by(AllExpData, Exp, WMSize, NSub) %>% 
  dplyr::summarize(m_WMCrr = mean(WMCrr), n =n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.

3 comparison between Exp.4a and Exp.4b)

3.1 correct rate

#plot WM correct rates
dplyr::group_by(AllExpData%>%filter(Exp %in% c('Exp4a', 'Exp4b')), Exp, WMSize, NSub, gap) %>% 
  dplyr::summarize(m_WMCrr = mean(WMCrr), n = n(), se_WMCrr = sd(WMCrr)/sqrt(n-1)) -> correctrate_gap
## `summarise()` has grouped output by 'Exp', 'WMSize', 'NSub'. You can override
## using the `.groups` argument.
write.csv(correctrate_gap, paste0(modelPath, '/rlt/correctrate_gap.csv'))

correctrate_gap%>%
  dplyr::group_by(Exp, WMSize, gap)%>%
  dplyr::summarize( n = n(),
             mean_WMCrr = mean(m_WMCrr), se_WMCrr = sd(m_WMCrr)/sqrt(n-1) ) -> meanForPlot_gap
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
WMCrr_gap <- ggplot(meanForPlot_gap, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr, group =interaction(Exp,gap), color = gap, shape = Exp, linetype = Exp))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  coord_cartesian(ylim = c(0.5, 1)) +
  colorSet5+
  labs(x = "Memory load", y = "Mean accuracy in WM task") +theme_new 
WMCrr_gap

ggsave(paste0(getwd(), "/figures/WMCrr_gap.png"), WMCrr_gap, width = 4, height = 4)

3.2 observed RP

AllExpData%>% filter(Exp %in% c('Exp4a', 'Exp4b')) %>%
  dplyr::group_by(Exp, curDur, WMSize, NSub, gap) %>% 
  dplyr::summarize(n = n(),
                   m_repDur = mean(repDur),
                   se_repDur = sd(repDur)/sqrt(n-1)) ->RP_obs_gap_bias
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize', 'NSub'. You can
## override using the `.groups` argument.
RP_obs_gap_bias$rep_bias = RP_obs_gap_bias$m_repDur - RP_obs_gap_bias$curDur


ezANOVA(data = RP_obs_gap_bias%>%filter(Exp =='Exp4b'), dv= rep_bias, wid=NSub, within= .(curDur, gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
## $ANOVA
##              Effect DFn DFd           F            p p<.05          ges
## 1            curDur   1  15 96.33440029 6.384821e-08     * 0.8176127533
## 2               gap   1  15  0.03549505 8.530873e-01       0.0001018033
## 3            WMSize   2  30  9.48732989 6.415750e-04     * 0.0895670700
## 4        curDur:gap   1  15 16.35270907 1.060464e-03     * 0.0220610259
## 5     curDur:WMSize   2  30  9.34453850 7.003993e-04     * 0.0280287517
## 6        gap:WMSize   2  30  2.51037940 9.816022e-02       0.0027254302
## 7 curDur:gap:WMSize   2  30  1.63535038 2.117823e-01       0.0021874657
RP_jasp <- RP_obs_gap_bias%>%filter(Exp =='Exp4b')%>% select(c( "curDur", "WMSize","NSub","gap", "rep_bias"))%>%  
   pivot_wider(names_from = c("curDur", "WMSize", "gap"), values_from = rep_bias, names_sep="")
## Adding missing grouping variables: `Exp`
#write.csv(RP_jasp, 'my_bias_rep_Exp5.csv')
RP_obs_gap  <- ggplot( data = RP_obs_gap_bias%>%
  dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur,  color=as.factor(WMSize), shape = gap)) +
  geom_point()+
  geom_line()+
  geom_abline(slope=1, intercept=0)+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+
   scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override
## using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_obs.png"), RP_obs_gap, width = 6, height = 4)

RP_obs_gap

3.3 central tendency effect (slope) and IPs (linear regression)

#Observed Indifference Point
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}

obs_InP <- AllExpData %>% 
  dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
obs_InP_Exp1_4 <- AllExpData %>% 
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
obs_InP_Exp1_4$model = NULL
obs_InP_Exp1_4$data = NULL
obs_InP_Exp1_4$inP = obs_InP_Exp1_4$Intercept /(1-obs_InP_Exp1_4$slope)

obs_InP_slope_Exp1_jasp <- obs_InP_Exp1_4 %>% filter(Exp == 'Exp1') %>% select(c("WMSize", "NSub","inP", "slope"))%>% pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp1_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp1_jasp.csv'))

obs_InP_slope_Exp2_jasp <- obs_InP_Exp1_4 %>% filter(Exp == 'Exp2') %>% select(c("WMSize", "NSub","inP", "slope"))%>% pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp2_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp2_jasp.csv'))

obs_InP_slope_Exp3_jasp <- obs_InP_Exp1_4 %>% filter(Exp == 'Exp3') %>% select(c("WMSize", "NSub","inP", "slope"))%>% pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp3_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp3_jasp.csv'))

obs_InP_slope_Exp4a_jasp <- obs_InP_Exp1_4 %>% filter(Exp =='Exp4a') %>% select(c("WMSize", "NSub", "inP","slope"))%>% 
  pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp4a_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp4a_jasp.csv'))

obs_InP_Exp4b_jasp <- obs_InP %>% filter(Exp == 'Exp4b') %>% select(c("WMSize", "NSub", "gap","inP", "slope"))%>% pivot_wider(names_from = c("WMSize", "gap"), values_from = c(inP,slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_Exp4b_jasp, paste0(modelPath, '/rlt/obs_InP_Exp4b_jasp.csv'))
mobs_InP = obs_InP %>% group_by(Exp, WMSize, gap)%>%
  dplyr::summarise(n=n(),
                   m_Intercept = mean(Intercept),
                   se_Intercept= sd(Intercept)/sqrt(n-1),
                   m_inP = mean(inP),
                   se_inP = sd(inP)/sqrt(n-1),
                   m_slope = mean(slope),
                   se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
plt_InP_linear_gap<- ggplot(data = mobs_InP, aes(x=WMSize, y=m_inP, group = gap, color = gap))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.3,  aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
  labs(colour = "Gap")+colorSet3+
  facet_wrap(~Exp)+
  xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
  theme(legend.position = "top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear_gap.png"), plt_InP_linear_gap, width = 5, height = 5)

plt_InP_linear_gap

ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= inP, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
##       Effect DFn DFd        F           p p<.05         ges
## 2        gap   1  15 4.065468 0.062041074       0.019543220
## 3     WMSize   2  30 8.762649 0.001006806     * 0.078253848
## 4 gap:WMSize   2  30 0.584382 0.563670241       0.003061788
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W           p p<.05
## 3     WMSize 0.4515805 0.003829536     *
## 4 gap:WMSize 0.5210587 0.010428128     *
## 
## $`Sphericity Corrections`
##       Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 3     WMSize 0.6458198 0.004991546         * 0.6808339 0.004255437         *
## 4 gap:WMSize 0.6761593 0.502620038           0.7194299 0.512209537
plt_RP_slope_linear_gap<- ggplot(data = mobs_InP, aes(x= WMSize, y=m_slope, group = gap,color = gap, shape = gap))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.3,  aes(ymin = m_slope - se_slope, ymax = m_slope + se_slope), position = position_dodge(width = 0.2)) +theme_new+
  facet_wrap(~Exp)+
  labs(colour = "Gap", shape = "Gap")+colorSet3+
  xlab(' ')+ylab("Slope of reproduction (s)")+
  theme(legend.position = "top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_RP_slope_linear_gap.png"), plt_RP_slope_linear_gap, width = 5, height = 5)

plt_RP_slope_linear_gap

ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= slope, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
##       Effect DFn DFd         F            p p<.05         ges
## 2        gap   1  15 16.352709 0.0010604639     * 0.027930881
## 3     WMSize   2  30  9.344539 0.0007003993     * 0.035428935
## 4 gap:WMSize   2  30  1.635350 0.2117822688       0.002784549
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W          p p<.05
## 3     WMSize 0.8669475 0.36808622      
## 4 gap:WMSize 0.6136631 0.03277255     *
## 
## $`Sphericity Corrections`
##       Effect       GGe       p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.8825716 0.001234014         * 0.9914082 0.0007299619         *
## 4 gap:WMSize 0.7213254 0.219115993           0.7775271 0.2179558985

4 load model results

4.1 check model parameters

### Average Parameters 
mm_Baypar <- dplyr::group_by(Bayparlist, Exp) %>%
  dplyr::summarize( m_sig_s2 = mean(sig_s2),
                    m_sig_pr2_log = mean(sig_pr2_log),
                    m_ks= mean(ks), 
                    m_kr = mean(kr), 
                    m_ls = mean(ls), 
                    m_ts = mean(ts),
                    m_mu_pr = mean(mu_pr),
                    m_sig_pr2 = mean(sig_pr2),
                    m_mu_pr_log= mean(mu_pr_log),
                    m_sig_mn2 = mean(sig_mn2),
                    n= n(),
                    se_sig_s2 = sd(sig_s2)/sqrt(n-1),
                    se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
                    se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
                    se_ks = sd(ks)/sqrt(n-1),
                    se_kr = sd(kr)/sqrt(n-1),
                    se_ls =sd(ls)/sqrt(n-1),
                    se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1)) 
mm_Baypar

4.2 Average Parameters

mBaypar_sub <- dplyr::group_by(Bayparlist, Exp, NSub) %>%
  dplyr::summarize( m_sig_s2 = mean(sig_s2),
                    m_sig_pr2_log = mean(sig_pr2_log),
                    m_ks= mean(ks), 
                    m_kr = mean(kr), 
                    m_ls = mean(ls), 
                    m_ts = mean(ts),
                    m_mu_pr = mean(mu_pr),
                    m_sig_pr2 = mean(sig_pr2),
                    m_mu_pr_log= mean(mu_pr_log),
                    n= n(),
                    se_sig_s2 = sd(sig_s2)/sqrt(n-1),
                    se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
                    se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
                    se_ks = sd(ks)/sqrt(n-1),
                    se_kr = sd(kr)/sqrt(n-1),
                    se_ls =sd(ls)/sqrt(n-1),
                    se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1)) 
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
fig_ks_subj = ggplot(mBaypar_sub, aes(Exp, m_ks, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  #facet_wrap(~Exp)+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor Ks')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_ks_subj.png"), fig_ks_subj, width = 6, height = 6)
fig_ks_subj

fig_kr_subj = ggplot(mBaypar_sub, aes(Exp, m_kr, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  #facet_wrap(~Exp)+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor Kr')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_kr_subj.png"), fig_kr_subj, width = 6, height = 6)
fig_kr_subj

fig_ts_subj = ggplot(mBaypar_sub, aes(Exp, m_ts, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor ts')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_ts_subj.png"), fig_ts_subj, width = 6, height = 6)
fig_ts_subj

fig_ls_subj = ggplot(mBaypar_sub, aes(Exp, m_ls, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  #facet_wrap(~Exp)+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor ls')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/mfig_ls_subj.png"), fig_ls_subj, width = 6, height = 6)
fig_ls_subj

fig_mu_pr_subj = ggplot(mBaypar_sub, aes(Exp, m_mu_pr, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'mean of prior')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_mu_pr_subj.png"), fig_mu_pr_subj, width = 6, height = 6)
fig_mu_pr_subj

fig_sig_pr2_subj = ggplot(mBaypar_sub, aes(Exp, m_sig_pr2, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'mean of variance')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_sig_pr2_subj.png"), fig_sig_pr2_subj, width = 6, height = 6)
fig_sig_pr2_subj

5 Prediction results

#load prediction
AllDat_predY <- read.csv(paste0(modelPath, "/rlt/AllDat_predY_",modelversion,".csv"))
AllDat_predY$WMSize <- as.factor(AllDat_predY$WMSize)
levels(AllDat_predY$WMSize) = c("low", "medium",  "high")
AllDat_predY$pred_Bias = AllDat_predY$mu_r - AllDat_predY$curDur
AllDat_predY$predErr = AllDat_predY$mu_r - AllDat_predY$repDur
AllDat_predY$relatErr = AllDat_predY$predErr / AllDat_predY$repDur
AllDat_predY[which(AllDat_predY$Exp == "Exp4"),"Exp"] = "Exp4a"
AllDat_predY[which(AllDat_predY$Exp == "Exp5"),"Exp"] = "Exp4b"
#calculate the mean reproduction biases for the five given intervals for all subjects
mpredY_sub <- dplyr::group_by(AllDat_predY, curDur, Exp, NSub, WMSize) %>%
  dplyr::summarize(n = n(), 
                   m_repDur = mean(repDur), 
                   sd_repDur = sd(repDur),
                   m_mu_r = mean(mu_r),
                   m_sig_r = mean(sig_r),
                   m_wp = mean(wp),
                   se_wp = sd(wp)/sqrt(n-1),
                   log_lik =mean(log_lik),
                   cv =sd_repDur/ m_repDur,
                   pred_cv = mean(sig_r/mu_r),
                   predRP_err = mean(m_mu_r-m_repDur),
                   predVar_err = mean(m_sig_r-sd_repDur),
                   predRP_rerr = mean(abs(m_mu_r-m_repDur)/m_repDur),
                   predVar_rerr = mean(abs(m_sig_r-sd_repDur)/sd_repDur),
                   predcv_err = pred_cv-cv,
                   predcv_rerr = mean(abs(pred_cv-cv)/cv))
## `summarise()` has grouped output by 'curDur', 'Exp', 'NSub'. You can override
## using the `.groups` argument.
mpredY_sub_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'm_repDur', 'curDur')), mpredY_sub$curDur) 
mpredY_sub_cv_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'cv', 'curDur')), mpredY_sub$curDur) 

mpredY_sub_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'curDur' ,'m_repDur')), mpredY_sub$WMSize) 
mpredY_sub_cv_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize','curDur', 'cv')), mpredY_sub$WMSize) 
mpredY_sub_jasp_RP = NULL
mpredY_sub_jasp_cv = NULL
mpredY_sub_jasp_RP1 = NULL
mpredY_sub_jasp_cv1 = NULL
for(i in 1: length(mpredY_sub_split)){
  temp = mpredY_sub_split[[i]]
  curDur = unique(temp$curDur)
  temp$curDur = NULL
  colnames(temp) = c('Exp', 'NSub', 'WMSize', paste0('m_repDur_', curDur))
  
  temp_cv = mpredY_sub_cv_split[[i]]
  curDur = unique(temp_cv$curDur)
  temp_cv$curDur = NULL
  colnames(temp_cv) = c('Exp', 'NSub', 'WMSize', paste0('cv_', curDur))
  if(i == 1){
    mpredY_sub_jasp_RP = temp
    mpredY_sub_jasp_cv = temp_cv
  }
  else{
    mpredY_sub_jasp_RP = left_join(mpredY_sub_jasp_RP, temp, by=c("Exp",  "NSub", "WMSize"))
    mpredY_sub_jasp_cv = left_join(mpredY_sub_jasp_cv, temp_cv, by=c("Exp",  "NSub", "WMSize"))
  }
  
}
for (i in 1: length(mpredY_sub_split1)){
  temp1 = mpredY_sub_split1[[i]]
  WMSize = unique(temp1$WMSize)
  temp1$WMSize = NULL
  colnames(temp1) = c('Exp', 'NSub', 'curDur', paste0('m_repDur_', WMSize))
  
  temp_cv1 = mpredY_sub_cv_split1[[i]]
  WMSize = unique(temp_cv1$WMSize)
  temp_cv1$WMSize = NULL
  colnames(temp_cv1) = c('Exp', 'NSub', 'curDur', paste0('cv_', WMSize))
  if(i == 1){
    mpredY_sub_jasp_RP1 = temp1
    mpredY_sub_jasp_cv1 = temp_cv1
  }
  else{
    mpredY_sub_jasp_RP1 = left_join(mpredY_sub_jasp_RP1, temp1, by=c("Exp",  "NSub", "curDur"))
    mpredY_sub_jasp_cv1 = left_join(mpredY_sub_jasp_cv1, temp_cv1, by=c("Exp",  "NSub", "curDur"))
  }
}

write_csv(mpredY_sub_jasp_RP, paste0(modelPath, '/rlt/RP_Bias_jasp.csv'))
write_csv(mpredY_sub_jasp_cv, paste0(modelPath, '/rlt/cv_jasp.csv'))
write_csv(mpredY_sub_jasp_RP1, paste0(modelPath, '/rlt/RP_Bias_jasp1.csv'))
write_csv(mpredY_sub_jasp_cv1, paste0(modelPath, '/rlt/cv_jasp1.csv'))
write_csv(dplyr::group_by(mpredY_sub, curDur, NSub) %>%
  dplyr::summarize(m_cv = mean(cv))%>%spread(curDur, m_cv), paste0(modelPath, '/rlt/m_cv.csv'))
## `summarise()` has grouped output by 'curDur'. You can override using the
## `.groups` argument.
mpredY_sub$RP_bias = mpredY_sub$m_repDur -mpredY_sub$curDur
mpredY_sub_new <- dplyr::group_by(mpredY_sub, curDur, Exp, NSub) %>%
  dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(curDur, m_RP_bias)
## `summarise()` has grouped output by 'curDur', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp1'), paste0(modelPath, '/rlt/RP_Bias_exp1.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp2'), paste0(modelPath, '/rlt/RP_Bias_exp2.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_exp3.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_exp4a.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_exp4b.csv'))


mpredY_sub_WMsize <- dplyr::group_by(mpredY_sub, WMSize, Exp, NSub) %>%
  dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(WMSize, m_RP_bias)
## `summarise()` has grouped output by 'WMSize', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp3.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4a.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4b.csv'))
#### predicted data
m_predY <- mpredY_sub%>% 
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_sd_repDur = mean(sd_repDur),
                   m_m_sig_r =mean(m_sig_r),
                   m_m_mu_r = mean(m_mu_r),
                   m_m_wp = mean(m_wp),
                   n = n(),
                   m_se_wp = sd(se_wp)/sqrt(n-1),
                   log_lik =mean(log_lik),
                   mpredRP_err = mean(predRP_err),
                   mpredVar_err = mean(predVar_err),
                   mpredRP_rerr = mean(predRP_rerr),
                   mpredVar_rerr = mean(predVar_rerr),
                   cv= mean(cv),
                   pred_cv = mean(pred_cv),
                   mpredcv_err = mean(predcv_err),
                   mpredcv_rerr = mean(predcv_rerr))
m_predY_acc =  mpredY_sub%>% 
  dplyr::group_by(Exp) %>% 
  dplyr::summarize(mpred_rerr = mean(predRP_rerr)*100,
                   mpredVar_rerr = mean(predVar_rerr)*100, 
                   mpredcv_rerr = mean(predcv_rerr)*100)
m_predY_acc

6 WAIC and LOO-CV

#extract waic and loo-cv from parameter list
m_WAIC <- dplyr::group_by(Bayparlist, Exp) %>%
  dplyr::summarize(n =n(),
                   m_looic = mean(looic),
                   m_waic = mean(waic),
                   se_waic = sd(waic)/sqrt(n-1),
                   se_looic = sd(looic)/sqrt(n-1),
                   m_p_loo = mean(p_loo),
                   m_elpd_loo = mean(elpd_loo),
                   m_se_looic = mean(se_looic),
                   m_se_p_loo = mean(se_p_loo),
                   m_p_waic = mean(p_waic),
                   m_se_waic = mean(se_waic)) 

m_WAIC
#load test results
AllDat_newY <- read.csv(paste0(modelPath, "/rlt/AllDat_newY_",modelversion,".csv"))
AllDat_newY$WMSize <- as.factor(AllDat_newY$WMSize)
levels(AllDat_newY$WMSize) = c("low", "medium",  "high")
AllDat_newY[which(AllDat_newY$Exp == "Exp4"), "Exp"] = "Exp4a"
AllDat_newY[which(AllDat_newY$Exp == "Exp5"), "Exp"] = "Exp4b"

7 Plot Results

7.1 RP biase

  RP_bias  <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY, aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp)) +
  labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
  colorSet3+guides(shape="none")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias.png"), RP_bias, width = 6, height = 6)

RP_bias

RP_bias_obs  <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%  
  dplyr::group_by(Exp, curDur, WMSize, NSub) %>% 
  dplyr::summarize(n = n(),
                   m_repDur = mean(repDur),
                   se_repDur = sd(repDur)/sqrt(n-1)) %>%
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur,  color=as.factor(WMSize))) +
  geom_point()+
  geom_line()+
  #geom_errorbar(width=.2,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+
   scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)

RP_bias_obs

RP_bias_obs  <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%  
  dplyr::group_by(Exp, curDur, WMSize, NSub) %>% 
  dplyr::summarize(n = n(),
                   m_repDur = mean(repDur),
                   se_repDur = sd(repDur)/sqrt(n-1)) %>%
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur,  color=as.factor(WMSize))) +
  geom_point()+
  geom_line()+
  geom_errorbar(width=.2,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+
   scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)

RP_bias_obs

RP  <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur,  color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY, aes(x=curDur, y=m_mu_r, color=WMSize)) +
  geom_abline(slope=1, intercept=0)+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+theme_new +theme(legend.position="top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP.png"), RP, width = 6, height = 6)

RP

7.2 RP Slope

# calculate the slope of the cv curve
RPslope_model <- function(df) {
  lm(m_repDur ~ log(curDur), data = df)
}

RPslopes <- mpredY_sub %>% 
  dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
  mutate(model = map(data, RPslope_model)) %>% 
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = `log(curDur)`)  
RPslopes$data <- NULL
RPslopes$model <- NULL
mRPslopes <- RPslopes%>% dplyr::group_by(WMSize, Exp) %>%
  dplyr::summarize(m_Intercept = mean(Intercept), 
                   m_slope = mean(slope),
                   n = n(), 
                   se_slope = sd(slope)/sqrt(n-1),
                   se_Intercept = sd(Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mRPslopes, aes(Exp, m_slope, ymin = m_slope - se_slope, ymax = m_slope + se_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+ 
  labs(x = "", y = TeX("Slope of RP"), color = 'Memory Load') +
  theme_new +theme(legend.position="top")

plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

7.2.1 Anova analysis on slopes of RP

RPslopes$WMSize = as.factor(RPslopes$WMSize)
ezANOVA(data = RPslopes, dv= slope, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd         F            p p<.05        ges
## 2        Exp   4  75  2.465821 5.210071e-02       0.11127986
## 3     WMSize   2 150 24.943347 4.446854e-10     * 0.01567449
## 4 Exp:WMSize   8 150  6.439500 3.455596e-07     * 0.01617814
## 
## $`Mauchly's Test for Sphericity`
##       Effect        W         p p<.05
## 3     WMSize 0.943444 0.1160104      
## 4 Exp:WMSize 0.943444 0.1160104      
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.9464714 1.164847e-09         * 0.9702459 7.594066e-10         *
## 4 Exp:WMSize 0.9464714 6.643447e-07         * 0.9702459 4.968822e-07         *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp1'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F         p p<.05         ges
## 2 WMSize   2  30 1.628257 0.2131416       0.002398051
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.8320858 0.2761702      
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05       HFe    p[HF] p[HF]<.05
## 2 WMSize 0.8562273 0.2172718           0.9557549 0.214484
ezANOVA(data = RPslopes %>% filter(Exp =='Exp2'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05        ges
## 2 WMSize   2  30 21.20007 1.822755e-06     * 0.07646424
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.8461525 0.3105563      
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.8666656 7.302883e-06         * 0.9698478 2.493639e-06         *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp3'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd         F         p p<.05        ges
## 2 WMSize   2  30 0.4137346 0.6648914       0.00152841
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.9800567 0.8684769      
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05      HFe     p[HF] p[HF]<.05
## 2 WMSize 0.9804466 0.6610089           1.126392 0.6648914
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4a'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 7.246942 0.002705914     * 0.01861911
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.9730381 0.8258648      
## 
## $`Sphericity Corrections`
##   Effect      GGe       p[GG] p[GG]<.05      HFe       p[HF] p[HF]<.05
## 2 WMSize 0.973746 0.002974256         * 1.117022 0.002705914         *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4b'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 8.551317 0.001151174     * 0.03802106
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.8930404 0.4529994      
## 
## $`Sphericity Corrections`
##   Effect       GGe       p[GG] p[GG]<.05      HFe       p[HF] p[HF]<.05
## 2 WMSize 0.9033753 0.001756187         * 1.019764 0.001151174         *

7.2.2 Anova analysis on Intercept of RP

ezANOVA(data = RPslopes, dv= Intercept, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd         F            p p<.05         ges
## 2        Exp   4  75 0.6157298 6.526476e-01       0.028876673
## 3     WMSize   2 150 4.6236001 1.125794e-02     * 0.005792574
## 4 Exp:WMSize   8 150 6.6943967 1.762289e-07     * 0.032641724
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W            p p<.05
## 3     WMSize 0.8020579 0.0002855036     *
## 4 Exp:WMSize 0.8020579 0.0002855036     *
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.8347649 1.634417e-02         * 0.8515168 1.573692e-02         *
## 4 Exp:WMSize 0.8347649 1.479632e-06         * 0.8515168 1.191853e-06         *

7.3 CV

curDurItem <- unique(m_predY$curDur)
RP_CV <- ggplot(data= m_predY, aes(x=curDur, y= cv, color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data = m_newY, aes(x=curDur, y= m_sig_r/m_mu_r, color=WMSize)) +
  facet_grid(~Exp) +
  labs(x="Duration (s)", y="CV", shape=" ", color = "Memory Load")+ theme_new+
  colorSet3+guides(shape="none")+theme(strip.text.x = element_blank())
RP_CV

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_CV.png"), RP_CV, width = 6, height = 6)

7.4 CV slope

# calculate the slope of the cv curve
cvSlope_model <- function(df) {
  lm(log(cv_obs) ~ log(curDur), data = df)
}

mpredY <- dplyr::group_by(AllDat_predY, curDur, WMSize, Exp, NSub) %>%
  dplyr::summarize(m_repDur = mean(repDur), 
                   sd_repDur = sd(repDur),
                   n = n(), 
                   sd_repDur = sd(repDur),
                   m_mu_r = mean(mu_r),
                   m_sig_r = mean(sig_r),
                   wp = mean(wp),
                   log_lik =mean(log_lik))
## `summarise()` has grouped output by 'curDur', 'WMSize', 'Exp'. You can override
## using the `.groups` argument.
mpredY$cv_obs <- mpredY$sd_repDur/mpredY$m_repDur

CVslopes <- mpredY %>% 
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%  # nested data
  mutate(model = map(data, cvSlope_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(obs_cv_Intercept = `(Intercept)`, obs_cv_slope = `log(curDur)`)  # rename columns
CVslopes$data = NULL
CVslopes$model = NULL
#change the table struction of slopes for spss
CVslopes_list <-split(CVslopes, CVslopes$WMSize) 
CVslopes_spss = NULL
for (i in 1: length(CVslopes_list)){
  temp = CVslopes_list[[i]]
  WMSize = unique(temp$WMSize)
  temp$WMSize = NULL
  colnames(temp) = c('Exp',  'NSub', paste0('obs_cv_Intercept_',WMSize), paste0('obs_cv_slope_',WMSize))
  if(i == 1)
    CVslopes_spss = temp
  else
    CVslopes_spss = left_join(CVslopes_spss, temp, by=c("Exp",  "NSub"))
}
mCVslopes <- CVslopes%>% dplyr::group_by(WMSize, Exp) %>%
  dplyr::summarize(n =n(),
                   m_obs_cv_Intercept = mean(obs_cv_Intercept), 
                   m_obs_cv_slope = mean(obs_cv_slope),
                   se_obs_CV_slope = sd(obs_cv_slope)/sqrt(n-1),
                   se_obs_CV_Intercept = sd(obs_cv_Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mCVslopes, aes(Exp, m_obs_cv_slope, ymin = m_obs_cv_slope - se_obs_CV_slope, ymax = m_obs_cv_slope + se_obs_CV_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+ 
  labs(x = "", y = TeX("Slope of CV"), color = 'Memory Load') +
  theme_new +theme(legend.position="top")

plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_CVslope.png"), plt_CVslope, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_CVIntecept<- ggplot(mCVslopes, aes(Exp, m_obs_cv_Intercept, ymin = m_obs_cv_Intercept - se_obs_CV_Intercept, ymax = m_obs_cv_Intercept + se_obs_CV_Intercept, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+ 
  labs(x = "", y = TeX("Intercept of CV"), color = 'Memory Load') +
  theme_new +theme(legend.position="top")

plt_CVIntecept
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Intercept.png"), plt_CVIntecept, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## Figures in the MS
fig_CV_slope_Intecept<-ggarrange(plt_CVslope, plt_CVIntecept, common.legend = TRUE, ncol=2, nrow=1,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_CV_slope_Intecept.png"), fig_CV_slope_Intecept, width = 8, height = 4)
fig_CV_slope_Intecept

## Figures in the MS
fig3<-ggarrange(RP_bias, RP_CV, common.legend = TRUE, ncol=1, nrow=2,  labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig3.png"), fig3, width = 6, height = 5)
fig3

8 Indifference Pointtemp_data = sumexpSub_diff %>%

8.1 Indifference Point(linear regression)

#Observed Indifference Point
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}

obs_InP <- AllExpData %>% #filter(curDur %in% c(1.1,1.4)) %>%
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
write.csv(obs_InP, file = paste0(modelPath, "/rlt/obs_InP.csv"))
mobs_InP = obs_InP %>% group_by(Exp, WMSize)%>%
  dplyr::summarise(n=n(),
                   m_Intercept = mean(Intercept),
                   se_Intercept= sd(Intercept)/sqrt(n-1),
                   m_inP = mean(inP),
                   se_inP = sd(inP)/sqrt(n-1),
                   m_slope = mean(slope),
                   se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_InP_linear<- ggplot(data = mobs_InP, aes(x= Exp, y=m_inP, color = WMSize))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.3,  aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
  labs(colour = "Memory Load")+colorSet3+
  xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
  theme(legend.position = "top")


ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear.png"), plt_InP_linear, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP_linear
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

RP_Exp2  <- ggplot(data = m_predY%>% filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur,  color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_segment(data =mobs_InP%>%filter(Exp == 'Exp2'), aes(x = 0.3, y = 0.3*m_slope+m_Intercept, xend = 2.2, yend = 2.2*m_slope+m_Intercept), arrow = NULL)+  ##dotted lines
  geom_point(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x = m_inP, y = 0.1, color=WMSize), shape=21) + ## Intersection 
  geom_segment(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x=m_inP, y = 0.1, xend = m_inP, yend = m_inP), arrow = NULL)+ ##vertical lines for Intersection
 #geom_line(data= m_newY%>% filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r, color=WMSize)) +
  geom_abline(slope=1, intercept=0, linetype = 2)+
  labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+ ylim(0,2)

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_Exp2.png"), RP_Exp2, width = 4, height = 4)

RP_Exp2

8.2 Indifference Point (bootstraps)

#### Estimated Indifference Point
mRep_model <- function(df) {
  lm(mu_r ~ curDur, data = df)
}

#Observed Indifference Point
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}

# Custom function to find predicted indifference point 
getPredInP_boot <- function(df, idx){
  vars <- c('NSub', 'Exp', 'WMSize')
  gp_vars = syms(vars) 
  slopes <- df[idx, ] %>% 
    dplyr::group_by(!!!gp_vars) %>% nest()  %>%  # nested data
    mutate(model = map(data, mRep_model)) %>%  # linear regression
    mutate(slope = map(model, broom::tidy)) %>%  # get estimates out
    unnest(slope, .drop = TRUE) %>% # remove raw data
    select(-std.error,-statistic, -p.value) %>%  # remove unnessary clumns
    spread(term, estimate) %>%   # spread stimates
    dplyr::rename(minRP = `(Intercept)`, slope = curDur)  # rename columns
  slopes$inP = slopes$minRP /(1-slopes$slope)
  return(slopes$inP)
}


# Custom function to find observed indifference point 
getRPInP_boot <- function(df, idx){
  vars <- c('NSub', 'Exp', 'WMSize')
  gp_vars = syms(vars) 
  slopes <- df[idx, ] %>% 
    dplyr::group_by(!!!gp_vars) %>% nest()  %>%  # nested data
    mutate(model = map(data, obs_model)) %>%  # linear regression
    mutate(slope = map(model, broom::tidy)) %>%  # get estimates out
    unnest(slope, .drop = TRUE) %>% # remove raw data
    select(-std.error,-statistic, -p.value) %>%  # remove unnessary clumns
    spread(term, estimate) %>%   # spread stimates
    dplyr::rename(minRP = `(Intercept)`, slope = curDur)  # rename columns
  slopes$inP = slopes$minRP /(1-slopes$slope)
  return(slopes$inP)
}
#calculate the bootstrapped 95% confidence intervals 
generateCI = FALSE # tag for generation CI
if(generateCI){
  cilist <- NULL
  for(expname in unique(AllDat_predY$Exp)){
    for(nsub in unique(AllDat_predY$NSub)){
      for(WMsize in unique(AllDat_predY$WMSize)){
        dat = AllDat_predY %>% filter(Exp == expname, NSub == nsub, WMSize ==WMsize)
        set.seed(100) 
        num = 1000
        bs_predInP <- boot(dat, getPredInP_boot, R = num) 
        bs_RPInP <- boot(dat, getRPInP_boot, R = num) 
        ci = data.frame(
          Exp = expname,
          NSub = nsub,
          WMSize = WMsize,
          sd_predInP_boot = sd(bs_predInP$t),
          m_predInP_boot = median(bs_predInP$t),
          sd_RPInP_boot = sd(bs_RPInP$t),
          mRP_InP_boot = median(bs_RPInP$t)
        )
        cilist = data.frame(rbind(cilist, ci))
      }
    }
  }
  write.csv(cilist, file = paste0("ci_list_median_1000.csv"))
}
# load the generated indifference point values and mark the outlier
cilist = read.csv(paste0("ci_list_median_1000.csv"))
cilist$Exp = as.factor(cilist$Exp)
cilist$WMSize= as.factor(cilist$WMSize)
cilist$inPOutlier = FALSE
cilist[which(cilist$mRP_InP_boot > 1.7 |cilist$mRP_InP_boot < 0.5 | cilist$m_predInP_boot < 0.5|cilist$m_predInP_boot > 1.7),"inPOutlier"] = TRUE

#check if the outlier is the same as the outliers in variable slope_pr
cilist%>% filter(inPOutlier == TRUE)

8.3 Inifference Point (curve)

AllDat_newY$predErr = AllDat_newY$mu_r - AllDat_newY$curDur

temp_newY <- AllDat_newY %>% filter(curDur > 0.8, curDur < 1.1) %>% select(Exp, WMSize, NSub, predErr, curDur)

InP_curve<- temp_newY%>% dplyr::group_by(Exp, WMSize, NSub)%>%
  dplyr::summarise(minErr = min(abs(predErr)), idx = which.min(abs(predErr)))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
InP_curve$InP_curve = temp_newY[InP_curve$idx,]$curDur 
InP_curve$y = temp_newY[InP_curve$idx,]$predErr + temp_newY[InP_curve$idx,]$curDur
InP_curve

8.4 Inifference Point (Equation)

#sig_sm2 =  sig_s2 + ls*(size[i]);
#wp_pr[i] = sig_sm2 /(sig_sm2 + sigma_pr2); 

Bayparlist$wp_1 = (Bayparlist$sig_s2 + Bayparlist$ls *1 )/ (Bayparlist$sig_s2 + Bayparlist$ls *1 + Bayparlist$sig_pr2_log)

Bayparlist$wp_3 = (Bayparlist$sig_s2 + Bayparlist$ls *2) / (Bayparlist$sig_s2 + Bayparlist$ls *2 + Bayparlist$sig_pr2_log)

Bayparlist$wp_5 = (Bayparlist$sig_s2 + Bayparlist$ls *3) / (Bayparlist$sig_s2 + Bayparlist$ls *3 + Bayparlist$sig_pr2_log)

#log(IP)=[kr*M + (1-w_p)ks*M+ w_p*mu_p]/w_p 

size = 1 
Bayparlist$InP_1 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_1)*Bayparlist$ks*size + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)
size = 2
Bayparlist$InP_3 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_3)*Bayparlist$ks*size + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)
size = 3
Bayparlist$InP_5 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_5)*Bayparlist$ks*size + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist$sig2_post_1 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5

Bayparlist$sig2_post_3 =(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5 

Bayparlist$sig2_post_5 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5

Bayparlist%>%filter(Exp == 'Exp2')%>%select("NSub", "Exp", "sig2_post_1", "sig2_post_3","sig2_post_5")
Bayparlist$InP_old_1 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_1)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5 + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)

Bayparlist$InP_old_3 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_3)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5 + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)

Bayparlist$InP_old_5 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_5)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5 + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist%>%select(NSub, Exp, wp_1, wp_3, wp_5, InP_1, InP_3, InP_5, InP_old_1, InP_old_3, InP_old_5)%>%
  unite(newcol1,  wp_1, InP_1, InP_old_1)%>%
  unite(newcol3,  wp_3, InP_3, InP_old_3)%>%
  unite(newcol5, wp_5, InP_5, InP_old_5)%>%
  melt(id.vars = c("NSub", "Exp")) %>%
  dplyr::rename(
    WMSize = variable,
    newcol = value
  )%>%
  separate(newcol,  c("wp", "InP_Eq", "InP_Eq_old"), sep = "_") -> Baypar_WMSize
Baypar_WMSize$WMSize <- substr(Baypar_WMSize$WMSize, 7, 7)
Baypar_WMSize$WMSize <- as.factor(Baypar_WMSize$WMSize)
Baypar_WMSize$InP_Eq = as.numeric(Baypar_WMSize$InP_Eq)
Baypar_WMSize$InP_Eq_old = as.numeric(Baypar_WMSize$InP_Eq_old)
levels(Baypar_WMSize$WMSize) = c("low", "medium",  "high")


##save csv for SPSS
left_join(Baypar_WMSize, InP_curve%>%select("Exp", "WMSize","NSub", "InP_curve"),  by = c("Exp","WMSize","NSub")) -> Baypar_WMSize

left_join(Baypar_WMSize, cilist %>% select("Exp", "NSub", "WMSize", "m_predInP_boot", "mRP_InP_boot"), by = c("Exp","WMSize","NSub"))-> Baypar_WMSize

8.5 check indifference Points

Baypar_WMSize_melt <- melt(Baypar_WMSize%>%select("Exp","WMSize","NSub", "InP_Eq", "InP_Eq_old", "InP_curve", "m_predInP_boot", "mRP_InP_boot"), id.vars = c("Exp","WMSize","NSub"),
                           variable.name = "Type",
                           value.name = "InP" )


Baypar_WMSize_melt %>%filter(Type== "InP_curve") %>% dplyr::group_by(Exp)%>%
  dplyr::summarise(n = n(),
                   m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))
Baypar_WMSize_melt$Exp = as.factor(Baypar_WMSize_melt$Exp)

Baypar_WMSize_melt%>% dplyr::group_by(Exp, WMSize, Type)%>%
  dplyr::summarise(n = n(),
                   m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))-> mBaypar_InP
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
#plot InP_curve(the intersections of the Prediction curve with the diagonal)
plt_InP<- ggplot(data = mBaypar_InP%>%filter(Type== "InP_curve"), aes(x= Exp, y=m_InP, color = WMSize, shape = Type))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  aes(ymin = m_InP - se_InP, ymax = m_InP + se_InP), position = position_dodge(width = 0.2)) +theme_new+
  labs(colour = "Memory Load")+colorSet3+
  xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
  theme(legend.position = "top")


ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP.png"), plt_InP, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

temp_inp = mBaypar_InP%>%filter(Type== "InP_curve", Exp == 'Exp2')
RP_bias_Exp2  <- ggplot(data = m_predY%>%filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY%>%filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  geom_point(data =temp_inp, aes(x = m_InP+0.03, y = -0.5, color=WMSize), shape=21) + ## Intersection 
  geom_segment(data =temp_inp, aes(x=m_InP+0.03, y = -0.5, xend = m_InP+0.03, yend = 0), arrow = NULL)+ ##vertical lines for Intersection
  labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
  colorSet3+guides(shape="none")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_Exp2.png"), RP_bias_Exp2, width = 5, height = 3)

RP_bias_Exp2

8.6 anova on observed InP

ezANOVA(data = Baypar_WMSize, dv= InP_curve, wid=NSub, within = .(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd        F            p p<.05        ges
## 2        Exp   4  75 1.122790 3.522214e-01       0.05010528
## 3     WMSize   2 150 4.862201 8.994962e-03     * 0.00766407
## 4 Exp:WMSize   8 150 5.529029 3.939989e-06     * 0.03393766
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W            p p<.05
## 3     WMSize 0.1284775 1.063161e-33     *
## 4 Exp:WMSize 0.1284775 1.063161e-33     *
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
## 3     WMSize 0.5343243 0.0280800589         * 0.535749 0.0279830198         *
## 4 Exp:WMSize 0.5343243 0.0004160714         * 0.535749 0.0004100909         *
pairwise.t.test(Baypar_WMSize$InP_curve, Baypar_WMSize$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Baypar_WMSize$InP_curve and Baypar_WMSize$Exp 
## 
##       Exp1  Exp2  Exp3  Exp4a
## Exp2  1.000 -     -     -    
## Exp3  1.000 1.000 -     -    
## Exp4a 1.000 1.000 1.000 -    
## Exp4b 0.125 0.030 0.049 0.690
## 
## P value adjustment method: holm
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp2'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 6.278925 0.005273799     * 0.04191182
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W            p p<.05
## 2 WMSize 0.1557357 2.221859e-06     *
## 
## $`Sphericity Corrections`
##   Effect       GGe      p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2 WMSize 0.5422216 0.02127963         * 0.5515802 0.02067639         *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp3'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 8.199485 0.001442758     * 0.04500834
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.03659746 8.793399e-11     *
## 
## $`Sphericity Corrections`
##   Effect       GGe      p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2 WMSize 0.5093199 0.01137853         * 0.5113321 0.01128146         *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4a'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F        p p<.05        ges
## 2 WMSize   2  30 1.507457 0.237775       0.01146743
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.09784573 8.586046e-08     *
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05       HFe     p[HF] p[HF]<.05
## 2 WMSize 0.5257197 0.2392211           0.5313463 0.2393698
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4b'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F          p p<.05        ges
## 2 WMSize   2  30 5.970484 0.00656511     * 0.08780338
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W            p p<.05
## 2 WMSize 0.1083146 1.749077e-07     *
## 
## $`Sphericity Corrections`
##   Effect       GGe      p[GG] p[GG]<.05       HFe     p[HF] p[HF]<.05
## 2 WMSize 0.5286291 0.02522631         * 0.5349067 0.0247749         *

9 Export data for spss

### change the table struction for spss
Bayparlist_list <-split(Baypar_WMSize %>% select(c("NSub","Exp","WMSize", "wp","InP_curve")), Baypar_WMSize$WMSize) 
Bayparlist_spss = NULL
for (i in 1: length(Bayparlist_list)){
  temp = Bayparlist_list[[i]]
  WMSize = unique(temp$WMSize)
  temp$WMSize = NULL
  colnames(temp) = c('NSub','Exp', paste0('wp_',WMSize), paste0('InP_curve_',WMSize))
  if(i == 1)
    Bayparlist_spss = temp
  else
    Bayparlist_spss = left_join(Bayparlist_spss, temp, by=c("Exp",  "NSub"))
}
write_csv(Bayparlist_spss, paste0(modelPath, '/rlt/Bayparlist_spss.csv'))

10 anova analysis

10.1 Anova on mean reproduction biases

mpredY_sub$RP_Bias = mpredY_sub$m_repDur-mpredY_sub$curDur
RP_bias_Anova <- ezANOVA(data = mpredY_sub, dv= RP_Bias, wid=NSub, within= .(curDur, WMSize), between = .(Exp), detailed = TRUE, return_aov = FALSE )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
RP_bias_Anova
## $ANOVA
##              Effect DFn DFd         SSn       SSd           F            p
## 1               Exp   4  75  0.20632235 6.7357598   0.5743293 6.820927e-01
## 2            curDur   1  75 47.70152785 7.0760974 505.5914872 4.589992e-35
## 4            WMSize   2 150  0.04568339 0.6961599   4.9216486 8.506722e-03
## 3        Exp:curDur   4  75  0.94184534 7.0760974   2.4956695 4.985449e-02
## 5        Exp:WMSize   8 150  0.25075426 0.6961599   6.7536818 1.507627e-07
## 6     curDur:WMSize   2 150  0.11975287 0.3666986  24.4927739 6.240642e-10
## 7 Exp:curDur:WMSize   8 150  0.12549340 0.3666986   6.4167178 3.670623e-07
##   p<.05         ges
## 1       0.013680912
## 2     * 0.762294526
## 4     * 0.003061808
## 3     * 0.059548049
## 5     * 0.016578279
## 6     * 0.007986470
## 7     * 0.008366110
# main effect of Duration  F(1.177, 3.532) = 377.965, p < .001, ηp² = .863.
(RP_bias_Anova$ANOVA)$DFn[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
#Duration × Experiment,  F(12, 240) = 2.506, p = .004, ηp² = .111
(RP_bias_Anova$ANOVA)$DFn[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
mpredY_sub <- ungroup(mpredY_sub)
res.aov <-  rstatix::anova_test(data = mpredY_sub, dv = RP_Bias, wid = NSub, within = c(curDur, WMSize), between = Exp)
## Warning: The 'wid' column contains duplicate ids across between-subjects
## variables. Automatic unique id will be created
get_anova_table(res.aov, correction = "GG")

10.2 variance of prior (anova)

ezANOVA(data = Bayparlist, dv= sig_pr2_log, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd         F       p p<.05        ges
## 1    Exp   4  75 0.5384571 0.70791       0.02791603
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd        SSn      SSd         F         p p<.05
## 1   4  75 0.04608107 1.121874 0.7701578 0.5480179

10.3 variance of motor noise (anova)

ezANOVA(data = Bayparlist, dv= sig_mn2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd        F         p p<.05        ges
## 1    Exp   4  75 0.602902 0.6617231       0.03115306
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd         SSn        SSd        F         p p<.05
## 1   4  75 0.001201703 0.03964056 0.568406 0.6863391

10.4 Ks anova

ezANOVA(data = Bayparlist, dv= ks, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 1    Exp   4  75 85.93971 3.114241e-27     * 0.8208993
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd       SSn       SSd        F            p p<.05
## 1   4  75 0.1584934 0.3319383 8.952723 5.765718e-06     *

10.5 mean of prior (anova)

ezANOVA(data = Bayparlist, dv= mu_pr, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd         F         p p<.05        ges
## 1    Exp   4  75 0.9353417 0.4482419       0.04751463
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd        SSn      SSd         F         p p<.05
## 1   4  75 0.02732829 1.139767 0.4495704 0.7723812

11 weight of prior

11.1 Weight of the prior \(w_p\)

Baypar_WMSize$wp = as.numeric(Baypar_WMSize$wp)
mwp <- Baypar_WMSize%>%dplyr::group_by(Exp, WMSize)%>% dplyr::summarise(m_wp = mean(wp), n= n(), se_wp = sd(wp)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_wp <- ggplot(mwp, aes(Exp, m_wp, ymin = m_wp - se_wp, ymax = m_wp + se_wp, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  #coord_cartesian(ylim = c(0.5, 1)) +
  colorSet5+
  labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load') +
  theme_new + theme(legend.position="top")

plt_wp
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_wp.png"), plt_wp, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

11.2 combine InP and wp

fig4<-ggarrange(plt_wp, plt_InP, common.legend = TRUE, ncol=2, nrow=1,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig4.png"), fig4, width = 6, height = 3)
fig4

fig5<-ggarrange(RP_bias_obs, plt_InP_linear, common.legend = TRUE, ncol=2, nrow=1, widths = c(5,3), labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig5.png"), fig5, width = 6, height = 3)
fig5

11.3 ANOVA analysis on wp

ezANOVA(data = Baypar_WMSize, dv= wp, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd          F            p p<.05        ges
## 2        Exp   4  75   6.864458 9.177058e-05     * 0.26663773
## 3     WMSize   2 150 189.342145 9.270731e-42     * 0.01709301
## 4 Exp:WMSize   8 150  33.257671 1.198007e-29     * 0.01207081
## 
## $`Mauchly's Test for Sphericity`
##       Effect          W            p p<.05
## 3     WMSize 0.02558321 1.242606e-59     *
## 4 Exp:WMSize 0.02558321 1.242606e-59     *
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.5064787 1.832257e-22         * 0.5067425 1.789171e-22         *
## 4 Exp:WMSize 0.5064787 4.021198e-16         * 0.5067425 3.954546e-16         *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp2'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 73.10886 2.92466e-12     * 0.03584474
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.03336632 4.604238e-11     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.5084831 3.075084e-07         * 0.5103134 2.944652e-07         *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4a'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 45.80731 7.62135e-10     * 0.02711716
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.01401012 1.059482e-13     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 2 WMSize 0.5035272 5.933032e-06         * 0.5042852 5.85177e-06         *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4b'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2 WMSize   2  30 72.16408 3.438001e-12     * 0.0546141
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.02522245 6.493984e-12     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.5063862 3.510759e-07         * 0.5077617 3.398985e-07         *
### pairwise.t.test on wp
pairwise.t.test(Baypar_WMSize$wp, Baypar_WMSize$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Baypar_WMSize$wp and Baypar_WMSize$Exp 
## 
##       Exp1    Exp2    Exp3 Exp4a
## Exp2  2.5e-13 -       -    -    
## Exp3  0.23    3.1e-08 -    -    
## Exp4a 0.23    4.1e-08 0.95 -    
## Exp4b 0.95    2.0e-11 0.69 0.69 
## 
## P value adjustment method: holm
##independent T test
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp1')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp1")))$wp
## t = 12.123, df = 47, p-value = 4.502e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2292634 0.3204897
## sample estimates:
## mean of x 
## 0.2748766
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp2')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp2")))$wp
## t = 15.978, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.5204341 0.6703593
## sample estimates:
## mean of x 
## 0.5953967
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp3')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp3")))$wp
## t = 17.739, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3138573 0.3941508
## sample estimates:
## mean of x 
## 0.3540041
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4a')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp4a")))$wp
## t = 11.932, df = 47, p-value = 7.948e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2968951 0.4173131
## sample estimates:
## mean of x 
## 0.3571041
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4b')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp4b")))$wp
## t = 11.53, df = 47, p-value = 2.664e-15
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2501705 0.3559257
## sample estimates:
## mean of x 
## 0.3030481

11.4 standard variance of Ds \(\sigma_{s}\)

ezANOVA(data = Bayparlist, dv= sig_s2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd        F          p p<.05       ges
## 1    Exp   4  75 3.409153 0.01287802     * 0.1538485
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd        SSn       SSd        F          p p<.05
## 1   4  75 0.01726871 0.1448521 2.235302 0.07314935
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(Bayparlist$sig_s2, Bayparlist$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Bayparlist$sig_s2 and Bayparlist$Exp 
## 
##       Exp1  Exp2  Exp3  Exp4a
## Exp2  0.078 -     -     -    
## Exp3  1.000 0.139 -     -    
## Exp4a 1.000 0.047 1.000 -    
## Exp4b 1.000 0.013 1.000 1.000
## 
## P value adjustment method: holm
##independent T test
t.test((Bayparlist%>%filter(Exp %in%c('Exp1')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp1")))$sig_s2
## t = 7.2054, df = 15, p-value = 3.046e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.02398402 0.04413435
## sample estimates:
##  mean of x 
## 0.03405919
t.test((Bayparlist%>%filter(Exp %in%c('Exp2')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp2")))$sig_s2
## t = 3.0831, df = 15, p-value = 0.007574
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.0247191 0.1354437
## sample estimates:
##  mean of x 
## 0.08008139
t.test((Bayparlist%>%filter(Exp %in%c('Exp3')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp3")))$sig_s2
## t = 11.682, df = 15, p-value = 6.232e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03170408 0.04585477
## sample estimates:
##  mean of x 
## 0.03877943
t.test((Bayparlist%>%filter(Exp %in%c('Exp4a')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp4a")))$sig_s2
## t = 6.5979, df = 15, p-value = 8.467e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.02042601 0.03992116
## sample estimates:
##  mean of x 
## 0.03017359
t.test((Bayparlist%>%filter(Exp %in%c('Exp4b')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp4b")))$sig_s2
## t = 4.5291, df = 15, p-value = 0.0003995
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.01172269 0.03256513
## sample estimates:
##  mean of x 
## 0.02214391
# calculate conhens d
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp1', 'Exp2')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp1", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.4080594
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp3')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.3966069
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp3', 'Exp4a')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp3", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.3980014
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp4a')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.4654873

11.5 cv slope

ezANOVA(data = CVslopes, dv= obs_cv_slope, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: Collapsing data to cell means. *IF* the requested effects are a subset
## of the full design, you must use the "within_full" argument, else results may be
## inaccurate.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd         F         p p<.05         ges
## 1    Exp   4  75 0.1126364 0.9777182       0.005971402
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd        SSn      SSd       F         p p<.05
## 1   4  75 0.08653032 1.378199 1.17722 0.3277357
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(CVslopes$obs_cv_slope, CVslopes$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  CVslopes$obs_cv_slope and CVslopes$Exp 
## 
##       Exp1 Exp2 Exp3 Exp4a
## Exp2  1    -    -    -    
## Exp3  1    1    -    -    
## Exp4a 1    1    1    -    
## Exp4b 1    1    1    1    
## 
## P value adjustment method: holm
##independent T test
t.test((CVslopes%>%filter(Exp %in%c('Exp1')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp1")))$obs_cv_slope
## t = -6.5917, df = 47, p-value = 3.406e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.2970793 -0.1581479
## sample estimates:
##  mean of x 
## -0.2276136
t.test((CVslopes%>%filter(Exp %in%c('Exp2')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp2")))$obs_cv_slope
## t = -7.47, df = 47, p-value = 1.59e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3369816 -0.1939864
## sample estimates:
## mean of x 
## -0.265484
t.test((CVslopes%>%filter(Exp %in%c('Exp3')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp3")))$obs_cv_slope
## t = -6.932, df = 47, p-value = 1.037e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3153826 -0.1735021
## sample estimates:
##  mean of x 
## -0.2444424
t.test((CVslopes%>%filter(Exp %in%c('Exp4a')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp4a")))$obs_cv_slope
## t = -5.4143, df = 47, p-value = 2.047e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3547829 -0.1625578
## sample estimates:
##  mean of x 
## -0.2586703
t.test((CVslopes%>%filter(Exp %in%c('Exp4b')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp4b")))$obs_cv_slope
## t = -7.4396, df = 47, p-value = 1.768e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3468257 -0.1991797
## sample estimates:
##  mean of x 
## -0.2730027
# calculate conhens d
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp1', 'Exp2')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.1178123
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp3')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.06429713
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp3', 'Exp4a')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.04045837
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp4a')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.01404826

12 Model prediction error

m_predErr_sub<- mpredY_sub%>% 
  dplyr::group_by(Exp, WMSize, NSub) %>% dplyr::summarise(
    mpredRP_err=mean(predRP_err),
    mpredVar_err=mean(predVar_err),
    mpredcv_err = mean(predcv_err),
    mpredRP_rerr = mean(predRP_rerr),
    mpredVar_rerr = mean(predVar_rerr),
    mpredcv_rerr = mean(predcv_rerr))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
m_predErr<- m_predY%>% 
  dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(
    mmpredcv_err = mean(mpredcv_err),
    mmpredRP_err=mean(mpredRP_err),
    mmpredVar_err=mean(mpredVar_err),
    mmpredRP_rerr = mean(mpredRP_rerr),
    mmpredVar_rerr = mean(mpredVar_rerr),
    mmpredcv_rerr = mean(mpredcv_rerr))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predErr
plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredcv_rerr*100, color = WMSize, alpha = .9)) + 
  #geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+ 
  geom_point() +
  geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredcv_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
  xlab('Relative prediction error of reproduction mean(%)')+ ylab('Relative prediction error of CV (%)')+colorSet3+
  facet_wrap(~Exp)+
  theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")

plt_rErrorScatter

plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredVar_rerr*100, color = WMSize, alpha = .9)) + 
  #geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+ 
  geom_point() +
  geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredVar_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
  xlab('Relative prediction error in the RP means (%)')+ ylab('Relative prediction error in the RP variance (%)')+colorSet3+
  facet_wrap(~Exp)+
  theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")

plt_rErrorScatter

plt_ErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_err, mpredVar_err, color = WMSize, alpha = .9)) + 
  geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+ 
  geom_point() +
  geom_point(data = m_predErr, aes(mmpredRP_err, mmpredVar_err, color = WMSize, alpha = .9, size = 1 ))+
  xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in the RP variance (s)')+colorSet3+
  facet_wrap(~Exp)+
  theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_ErrorScatter.png"), plt_ErrorScatter, width = 7, height = 4)

plt_ErrorScatter

13 model comparison (logarithmic vs. linear)

m_predErr_sub$model = 'logarithmic'
m_predErr$model = 'logarithmic'
linear_model = 'linear_rstan'
m_predErr_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_", linear_model, ".csv"))
m_predErr_linear$X = NULL
m_predErr_sub_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_sub_", linear_model, ".csv"))
m_predErr_sub_linear$X = NULL

m_predErr_sub_all = rbind(m_predErr_sub, m_predErr_sub_linear) 
m_predErr_all = rbind(m_predErr, m_predErr_linear)
m_predErr_all$WMSize = as.factor(m_predErr_all$WMSize)
levels(m_predErr_all$WMSize) =  c("low", "medium",  "high")
temp = m_predErr_all %>% filter(model == 'logarithmic') %>%summarise(abs_mmpredcv_err = abs(mmpredcv_err)) 

plt_Err_CV_all = ggplot(m_predErr_all, aes(abs(mmpredRP_err), abs(mmpredcv_err), group = interaction(model, Exp, WMSize), color = model, shape = Exp, size = WMSize)) + 
  geom_point() +
  geom_hline(yintercept = round(max(temp$abs_mmpredcv_err), 4)+0.0005, linetype='dashed')+
  xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in CV')+colorSet3+
  scale_shape_manual(values = c(6, 7, 16, 17,8)) +
  theme_new+ 
  theme(legend.position = 'top')+
 labs(size = 'Memory Load')+
  guides(colour = guide_legend(order = 1, nrow=2,byrow=TRUE),
         shape = guide_legend(order =2, nrow=2,byrow=TRUE),
            size = guide_legend(order = 3, nrow=3,byrow=TRUE))

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Err_CV_all.png"), plt_Err_CV_all, width = 7, height = 4)
## Warning: Using size for a discrete variable is not advised.
plt_Err_CV_all
## Warning: Using size for a discrete variable is not advised.

m_predY_acc =  m_predErr_sub_all%>% 
  dplyr::group_by(Exp, model) %>% 
  dplyr::summarize(mmpredRP_rerr = mean(mpredRP_rerr)*100,
                   mmpredVar_rerr = mean(mpredVar_rerr)*100,
                   mmpredcv_rerr = mean(mpredcv_rerr)*100,
                   mmpredRP_acc = (1-mean(mpredRP_rerr))*100,
                   mmpredVar_acc = (1-mean(mpredVar_rerr))*100,
                   mmpredCV_acc = (1-mean(mpredcv_rerr))*100)
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predY_acc